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Summary of Work Done Under Grant Sponsorship 

The Tnain purpose of this research has been to develop a rigorous theory and corre- 
sponding computational algorithms for through-the-thickness analysis of composite plates. 
This type of analysis is needed in order to find the elastic stiffness constants for a plate and 
to post-process the resulting plate solution in order to find approximate three-dimensional 
displacement, strain, and stress distributions throughout the plate. This also requires the 
development of finite deformation plate equations which are compatible with the through- 
the-thickness analyses. 

After about one year’s work, we settled on the variational-asymptotical method 1 
(VAM) as a suitable framework in which to solve these types of problems. VAM was 
applied to laminated plates with constant thickness in the work of Atilgan and Hodges 2 . 
The corresponding geometrically nonlinear global deformation analysis of plates was de- 
veloped by Hodges, Atilgan, and Danielson 3 . A different application of VAM, along with 
numerical results, was obtained by Hodges, Lee, and Atilgan^. An expanded version of 
this last paper 5 Has been submitted for publication in the AIA A Journal (Copies of these 
papers have been delivered to Mr. Hinnant.) One last paper was just completed and a 
copy is included as an appendix to this report. Summaries of the progress in the various 
categories we worked on follow. Technical details are in the papers. 


Work on Finite Deformation of Plates 

In Ref. 3, a set of kinematical and intrinsic equilibrium equations are derived for large 
deflection anH rotation but with small strain. This work showed that the drilling type 
rotation is not an independent variable in plate theory. If it is to be included in plate 
equations at alL there must be a constraint enforced between it and its definition in terms 
of other plate variables. Unless this constraint is enforced, the theory including the rotation 
as a separate variable is not valid. The main contribution of this paper is a complete set 
of geometrically exact strain-displacement relations for large deformation of plates. 

Also, the relationship between the drilling rotation and the other kinematical variables 
gives new in<rigHt into the dr illing moment and its role in beam-plate connectivity. An 
applied drilling moment at a point on a plate is not resisted at all by the plate. Such a 
moment, in order to have any physical resistance from the plate, must be applied over a 
finite area. Other than this, a point drilling moment can only be resisted by a plate if 
the plate model is derived from couple-stress elasticity 7 . Further study on this problem 
was shelved, because of the need to interpret the drilling rotational degrees of freedom in 
certain plate znd shell finite element formulations. There is an apparent conflict between 
the theoretical result of a plate having zero stiffness in response to a point drilling moment 
on one hand, and the finite element results found in the literature for the response of 
plates to point drilling moments. The papers in which these results appear are of several 
researchers, some of whom are highly regarded by the international community and by 
the principal investigator. A dialog has been initiated on this subject with some of these 
people, but to date we have no resolution to this question. 
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Work on Modeling Laminated Plates 

By modeling here we mean the calculation of elastic constants for a two-dimensional 
model based on known material constants for the three-dimensional body. This can only 
be done approximately, of course, and asymptotical methods are natural. In Ref. 2, the 
“first approximation” is exactly the same as classical laminated plate theory. The “sec- 
ond approximation” takes transverse shear deformation into account and is asymptotically 
correct only for plates with certain restrictions in their construction. To remove the re- 
strictions one must “kill” certain interaction terms in the strain energy, but the means for 
doing so for general laminated plates were not given in that paper. Indeed, the means for 
doing this were unknown to us at that tune. 

The development in Refs. 4 and 5 includes transverse shear in the “first approxima- 
tion” and is stopped there. Results from this theory were compared, for the cylindrical 
bending cas e, with results from the exact solution of Pagano for both cross-ply and shear- 
coupled 9 laminat ed plates. The resulting theory, termed a “neo-classical theory (NCPT), 
is at least as good as classical plate theory (CPT) in every case; for most cases NCPT 
is superior to CPT. For example, NCPT does a much better job on calculation of plate 
(two-dimensional) displacement than CPT. Also, many three-dimensional quantities are 
calculated much more accurately than could be achieved with CPT. Results for several 
different types of plates may be found in Ref. 5. 

As pointed out in Ref. 5, there is a need to develop higher approximations in order to 
(1) improve the overal performance of the theory for applications to thick plates and (2) 
improve the accuracy of predicted transverse strains and stresses, especially in situations 
when integration of the three-dimensional equilibrium equations cannot be accomplished, 
such as in the geometrically nonlinear case. 


Work on a Higher Approximation for NCPT 

In order to improve the asymptotical accuracy of NCPT one needs to introduce “de- 
grees of freedom” of the normal line element which will eliminate or “kill” interaction 
terms of the type identified in Ref. 2. In our work during the last reporting period we 
applied Sutyrin’s eigen-principle (described in the last report) to develop a refined theory 
for laminated plates of constant thickness. This theory includes CPT as a special case 
and even in its most elementary extension should surpass NCPT in predictive capability. 
Furthermore, because it is based on a variable number of “degrees of freedom” for the 
normal lin e element, it should allow users of future finite elements based on this theory to 
decide which of these degrees of freedom they want in their models. In principle, one can 
approach the accuracy of three-dimensional elasticity to any degree desired. The theory is 
completed as far as the derivation of the stiffness model based on eigenfunctions associated 
with a certain Stunn-Liouville operator. A computer code has been written and validated 
for the calculation of these eigenfunctions. Another code is presently under development 
which will use these eigenfunctions to build the stiffness model. For the details of the 
stiffness model, see Ref. 6. 
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Depending on what funding we can find, in the future we intend to generate results 
from this new theory for the cylindrical bending problem. This will allow us to inexpen- 
sively evaluate the number and types of normal line element degrees of freedom required for 
accurate solutions. The last step would then be to build a prototype plate finite element. 
We believe that this element would be far superior to any extant plate finite elements. 
We would want to work with someone who already has experience in two-dimensional fi- 
nite elements. The strength of our element would be in the theoretical foundation, not 
in the finite element technology itself. Dr. Sutyrin also has devised a means by which a 
discretized model based on three-dimensional finite elements can be reduced directly to 
plate elements. While this should be equivalent to our present semi-analytical approach, 
it may prove to be more convenient in certain contexts. 
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A Variable-Order Laminated Plate Theory 
Based on the Variational-Asymptotical Method 

Bok W. Lee; Vladislav G. Sutyrin;* and Dewey H. Hodges* 
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Abstract 

The variatfotiahasymptotical method is a mathemat- 
ical technique by which the three-dimensional analysis of 
laminated plate deformation can be split into a linear, one- 
dlm ensiooal, through-tbe-thkkness analysis and a nonlin- 
ear, two- dim ensional, plate analysis. The elastic constants 
used in tl* plate analysis are obtained from the through- 
t he- thickness analysis, along with approximate, closed-form 
three-dirnafifocal distributions of displacement, strain, and 
stress. In this paper, a theory based on this technique 
is developed which is capable of approximating three- 
dimensional elasticity to any accuracy desired. This theory 
is not developed using any of the usual approaches of lam- 
inated plate theory. That is, it is not based on any power 
series expansion through the thickness, nor is it based on 
introduction cf a set of variables which describe displace- 
ment in separate layers of laminated plates. Rather, the 
asymptotical method allows for the approximation of the 
through-d^-thickness behavior in terms of the eigenfunc- 
tions of a cstain S t urm-Lkxrville problem associated with 
the thickness coordinate. These eigenfunctions contain all 
the necessary information about the nonhomogeneities along 
the thicknes coordinate of the plate and thus possess the ap- 
propriate discontinuities in the derivatives of displacement. 
The theory is presented in this paper along with numerical 
results for the eigenfunctions of various laminated plates. 

Introduction 

For aerospace structures, laminated composite materi- 
als provide excellent opportunities for structural simplicity 
as well as elastic couplings for potential optimization of de- 
sign criteria. Although plates made of such materials have 
been used for some time in a variety of engineering appli- 
cations, simple and efficient methods for analyzing plates 
with anis otropy and nonhomogeneity are still needed. This 
is partly doe to rapid changes taking place in manufactur- 
ing te chnofo gy for composite materials and partly to ever- 
increasing demands for accuracy and efficiency. Much of 
what is done is based on classical plate theory (CPT) which, 
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although adequate for many engineering applications, has 
wefl known limitations due to the Kirchhoff hypothesis. 

Background 

M»ny attempts have been made to improve classical 
theor y by taking mtn account non-classical behavior such as 
transverse shear deformation and transverse normal stresses. 
From the time laminated fiber-reinforced composites were 
first introduced, numerous works have been published, the 
objectives of which include the improvement of CPT for lam- 
inated plate applications. This subject is discussed at length 
in review papers. 1 * 2 There are two main classes of methods 
for improving plate theory found in the literature: (1) P ower 
Series Methods: expansion of the displacement field variables 
into higher-order power series in the thickness coordinate; 
and (2) Layerwise Variables Methods: incorporation of sep- 
arate sets of displacement variables for each layer. Both of 
thpgp methods have known shortcoming. For example, no 
power series expansion can possibly render accurate results 
for quantities which may possess discontinuities, such as cer- 
tain components of strain and stress in laminated plates. 
The layerwise variables methods rely on a significant in- 
crease in the number of unknowns, a number which depends 
on the number of layers in the plate. A third method has 
received some attention in recent years, which involves an 
assumed displacement field with discontinuities allowed in 
through-the- thickness derivatives?* 4,5 There is no question 
that this method yields excellent results in some cases, but it 
jarW a systematic basis for choosing the displacement func- 
tions, and it does not yield an asymptotically correct result 
in general. 

Ref. 6 undertook a quite different approach. It does not 
involve a power-series expansion through the plate thick- 
ness, nor does it involve layerwise unknowns. Rather, the 
three-dimensional energy of a laminated plate was approx- 
imated following Berdichevsky’s variational- asymptotical 
methodology? Normally asymptotical methods are employed 
for analytical developments, but here such a method was 
used in a sort of semi-analytical approach. Namely, the 
theory leads to a Reissner-like plate theory, along with a 
set of elastic constants; it also provides a set of influence 
functions from which approximate three-dimensional dis- 
placement, strain, and stress fields can be determined once 
the plate equations are solved. The plate equations can be 
solved by any method desired, such as a two-dimensional 
finite element method. The analysis was restricted to lami- 
nated plates for which each lamina exhibits monoclinic m ar- 
terial symmetry about its middle surface. Their first ap- 
proximation is asymptotically correct for this case and coin- 
cides with classical laminated plate theory. However, their 
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second approximation is asymptotically correct only when 
each element of tbe reduced stiffness matrix Q (see Jones 8 ) 
is constant through the thickness of the entire plate. Al- 
though the theory is not asymptotically correct otherwise, it 
was intended f ? application to laminated composite plates. 
The limitation of their work stems from a term in the strain 
energy which was neglected. A suitable method to make 
this term vanish rigorously is needed in order to make the 
theory asymptotically correct, but was not given. 

A similar, but somewhat improved, approach was un- 
dertaken by Refs. 9 and 10, in which the estimation pro- 
cedure of R tL 6 was slightly modified to include transverse 
shear terms in the first approximation. Plates with cross- 
ply stacking sequences under cylindrical bending were taken 
as example problems in Ref. 9. In a later extension of this 
work? 0 plates with arbitrary stacking sequences undergoing 
cylindrical bending were taken as example problems. Tbe 
material configurations of these latter plates are not as sim- 
ple as of bi-directional plates, because of the influence 
of the coupling of transverse shear terms. Tbe distributions 
of three-dimensional displacement, strain, and stress were 
investigated for both cases by comparing the corresponding 
three-dimensional exact elasticity solution? 1,12 The theory 
of Refs. 9 and 10, termed tbe “neo-classical” plate theory 
(NCPT), was shown to be more accurate than CPT for thick, 
laminated plates; also, it was shown to yield results which 
are somewhat better than those of the theory of Ref. 6. 

Stiff, there were results reported for which the corre- 
lation of NCPT with the exact solution is not good. For 
example, when shear coupling exists, NCPT shows signif- 
icantly better correlation with the exact solution than for 
the bi-directional cases. It is necessary, then, to extend the 
validity of the theory to a higher approximation. Although 
it is not discussed in Refs. 9 and 10, such an extension re- 
quires certain interaction terms vanish. These terms 
are analogous to the one neglected in Ref. 6. For this reason 
a method for generalization of the theory of Refs. 9 and 10 
has been developed, and that is the subject of the present 
paper. 

Present Approach 

Tbe essence of the new approach, which guarantees dis- 
appearance of the interaction terms discussed above, is the 
introduction of new “degrees of freedom” into the three- 
dimensional displacement field. We are not using the term 
“degrees of freedom” in its usual sense. Here we mean plate 
(t.e., two-dimensional) displacement variables which are as- 
sociated with a particular deformation mode of the normal 
line element through the thickness. The shape functions 
of these new degrees of freedom are chosen as the eigen- 
functions associated with a certain Sturm-Liouville problem 
formulated by Sutyrin. 13 By choosing the associated warp- 
ing displacement to be orthogonal to the shape functions for 
each of tbe new degrees of freedom, the displacement field is 
uniquel v defined. The choice of shape functions from these 


eigenfunctions guarantees that any additional warping in- 
duced by the new degrees of freedom is of a higher order rel- 
ative to the new degrees of freedom themselves. The order of 
these new degrees of freedom relative to the strain depends 
cm the loading and the material constants. After obtaining 
the eigenvectors by using a one-dimensional finite element 
method, the plate elastic constants can be obtained through 
the variational-asymptotical method. Utilizing global defor- 
mation equations along with this resulting plate constitutive 
law, we will complete the formulation of the theory. 

In thin paper we first provide the details for tbe kine- 
matics of plate deformation in terms of classical plate dis- 
placement variables and the new degrees of freedom. Tbe 
strain field needed to develop a geometrically nonlinear plate 
theory is written in terms of these displacement variables. 
The small parameters are identified as f , the maximum 
strain in the plate, and 7, where h is the plate thick- 
ness and l is the characteristic length over which the de- 
formation varies in the deformed plate. The variational- 
asymptotical method is then used, along with a Ritz-type 
approximation of three-dimensional displacement variables 
in the through- the- thickness coordinate, to approximate tbe 
three-dimensional strain energy of the plate with a function 
of two-dimensional quantities only. The above-mentioned 
Sturm-Liouville problem is identified, the eigenfunctions of 
which guarantee the warping to be of higher order than 
any retained degree of freedom if they are chosen as the 
shape functions for the retained degrees of freedom. The 
two-dimensional strain energy function is then given as a 
function of material constants and eigenfunctions. 

We are able to calculate the exact solution for the 
Sturm-Liouville problem only for one- and two-layer plates. 
Thus, it was necessary to develop an approximate method 
of solution so that the plate theory could be finalized. We 
present a finite element analysis in one dimension (through 
the thickness) which we use to solve the Sturm-Liouville 
problem based on the shape functions of Ref. 14. Eigenval- 
ues and eigenfunctions for various laminated plates obtained 
by a one-dimensional finite element method are presented 
and, when possible, compared with the exact solution. Ap- 
plication of the theory will be addressed in a later paper. 

Theoretical Development 

The objective is to derive a strain energy function of a 
plate in terms of two-dimensional quantities only. It can be 
done only if some small parameters are present. We suppose 
that the parameters mentioned above, € and 7, are small. 

Three-Dimensional Formulation 

To begin we will formulate a three-dimensional devel- 
opment which shall be considered the exact solution to the 
plate problem. 

Undeformed State of Plate A typical point in the unde- 
formed plate can be located by introducing a Cartesian coor- 
dinate system X* in such a way that x a : {xi , x?) = x denotes 
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lengths ai-:cg orthogonal straight lines in the mid-surface of 
the un derrated plate, and X 3 = y = is the distance in 
the normal direction, where — 5 < s ^ 5 * Throughout the 
analvsis. Greek indices assume values 1 or 2; Latin indices 
assume values 1 , 2, and 3; and repeated indices are summed 
over their ranges. 

The spatial position vector f(xi,i 2 ,y) to 80 arbitrary 
point in t be undeformed plate can be written as 

r(x, y) = r(x) + yb(x) ( 1 ) 

where r (x) ts the spatial position vector of points on the 
mid-surface of a plate and b (x) is the unit normal vector. 
We will also need notation for unit rectors b a s ^ which, 
together with b, form an orthonormal triad. 

Since the variable y (and () is chosen specifically so that 
the spatial position vector r to a point on the reference mid- 
surface is the average position of points along the normal 
line at a particular value of Xi and xj. then 

r(x) = (r(x.O) ( 2 ) 

where the notation 



is used throughout the paper. 

Deformed State of Plate Without any restrictions the posi- 
tion rector X 2 , y) to an arbitrary point in the deformed 

plate can be represented by 

R(x,y) = R(x) + yB(x) + t^(x, C)B„(x) (4) 

where R(x) is the position vector of points on the deformed 
reference surface, B„(x): (Bi(x), B^x), Ba(x) = B(x)} 
is the refee nee orthonormal triad with vector B being or- 
thogonal to the deformed reference surface, and u n (x, C) are 
componffi-s of the general warping displacement of an arbi- 
trary' pocit in the deformed normal line, consisting of both 
in- and ccn-of-plane components so that all possible defor- 
mations are considered. 

The warping v n could not be defined uniquely as a func- 
tion of £ with an arbitrary choice of R(x) unless they are 
subject tc the constraints 

<v n (*,0)=0 (5) 

which means that 

R(*)=(r<x,C)) (6 ) 

The orthogonality of vector B to the reference surface 
means 

R a B = 0 (7) 


where the notation ( ) t « denotes the partial derivative with 
respect to x a . 


In order to eliminate the arbitrary rotation of vectors 
B Q around normal B we impose the following constraint 

B x • R . 2 = B 2 R.i ( 8 ) 


A schematic of the plate deformation is shown in Fig. 1. 

Uxkformod State Defocmcd State 



Fig. 1: Schematic of plate deformation 
(u is the displacement of the reference surface) 


Thus, Eq. (4) provides a convenient way to represent 
the arbitrary function R(x,y). The orientation of the triad 
Bn is now specified uniquely. 


Strain Field As shown in Ref. 6 , under the condition of 
small local rotation, the Jaumann strain components can be 
arranged in a 3 x 3 symmetric matrix T * , given by 


r* = i(x + x T ) - -f 

Xo =B, - j* 


(9) 


where / is the 3 x 3 identity matrix. 

Substituting Eq. (4) into the Jaumann strain, Eq. (9), 
one can express the strain field as a 6 x 1 column matrix 

r = in, 2 H, 2r ; 3 21^ r^j T so that 

r = |r h v' + r,7 + r/v I (10) 

n 
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where ^ and matrices Th , T t and I* are 

r< = 


r* = 


ro o o 
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0 o o 

1 o o 
o 1 o 

lo o 1J 


i <i 

0 0 


T t = [it r 2 ] 
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1 

0 

0 
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0 
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0 
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.0 

0 
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.0 

0 

0. 


’-{«} <-{x} (i2) 

Here we wiB introduce the twx>-dimensionaJ strain measures 
for in plane stretching and shear 


( R 4 B x -1 ) 

< R t i - Ba 4- R 4 * Bi > 

{ R 4 B*-l J 


(13) 


and the two-dimensional strain measures for bending and 
twist 


f B 4 -B 1 ) 

K = { B ! B 2 4 B a B x \ 

{ B 4 -Ba J 

The <gnaH parameter s can be now specified as 
E — max | 7 | 


(14) 


(15) 


A few rwnlinear terms in the strain field, which couple v 
and 7 , hav^ been neglected in Eq. ( 10 ) because the physically 
linear plate theory is considered. The form of the strain field 
is of great importance because h is now linear in 7 , t/, and 

This is the only point where £ as a small parameter is 
taken into account. 


Strain Erjzrgy of Plate The strain energy density per unit 
area for a plate can be written as 

(7 = i(r r Dr) (16) 

where D is the 6 x 6 symmetric material matrix in the B. 
basis. 

The three-dimensional Jaumann stress Z y which is con- 
jugate to the Jaumann strain T is 


Basic Three-Dimensional Problem 

The basic three-dimensional problem can now be rep- 
resented as the following minimization problem 


/Cl( 7 (x),f'(r,(),t) iI (i,C)) hdxidx 2 + 
4 terms with external forces — ► min 


(18) 


where the minimum should be found with respect to three 
arbitrary functions R(x), through which t(x) is calculated, 
and to three functions v„(i, £) which are subject to the con- 
straints of Eq. (5). 

Note that the orthonormal triad B n is not an inde- 
pendent variable, since it is & subject to the constraints in 
Ee& ( 7 ) and ( 8 ). It is completely determined by the function 
R(x). 


Dimensional Reduction 


Splitting of the Problem Now the functional space of all the 
admitted functions R(x, y) is split into sublayers with a 
choice of the x-dependent functions R(x). In each layer 
the functions tz n (x, £) are arbitrary under the constraint in 
Eq. (5). 

We ran solve this problem in two steps. Fust we are 
going to find functions u n (x,<) for any prescribed choice 
of functions R(x). As a result, we will have u„(x,C) as a 
functional of R(x) and C, and the functional, Eq. (18), will 
become dependent only on R(x). That functional will give 
us a two-dimensional plate theory. The second step is to 
solve that theory. 

Since the energy density U depends not only on func- 
tions v n (x, C) hut also on their derivatives with respect to x, 
then the result of the first step will be very complicated (it 
will contain a non-local dependence on R(x) in the general 
case) and cannot be obtained in an appropriate form unless 
we take advantage of some small parameters. 


Small Parameters Let us consider the situation where pa- 
rameters h, £ and e are present. Since no coefficient matrix 
of T, Eq. (11), depends on these parameters, it is clear that 
the first term of the expression for strain, Eq. (10), has or- 
der the second term has order £, and the third one has 
order The third term has order £ times that of the 
first. We should neglect this as a higher order term in the 
first approximation if we are going to expand the solution 
with respect to the small parameter In the future, this 
important circumstance allows us to avoid the presence of 
derivatives of unknown functions v n (x, £) with respect to x 
in the problem for any approximation and, hence, to solve 
it in an appropriate form. 


2 ? — jjy (17) Since our ttiaiti problem has become linear with respect 

to the unknown functions v n (x, £) and the two-dimensional 
strain measure 7 , the smallness of parameter e does not 
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need to be considered any more. This fact has already taken 
into account (see above under Strain Field ). We will expand 
the warping as a series with respect to the small 

parameter j only. In the absence of other small parameters, 
expansion in * is the same as expansion in h. That is why 
we can consider h to be the small parameter in spite of its 
dimension. 


D iscret ixation 

The problem may be solved numerically by discretizing 
with respect to the thickness variable C- Now the unknown 
functions are represented as the product of a matrix 

of shape functions 5(C) and a column matrix of nodal values 
of v(x, C)» which we will denote V 

v(x,C) = S(QV(x) (19) 


Substituting the discretized unknown function in 
Eq. (19) into the energy density, Eq. (16), while taking into 
account the strain, Eq. (10), one obtains 


2U = (tf\ rT EV+ (£) 2 V 7 [D U 7 + D hx V x ] + (2Q) 

+(l)(r r Z>«7 + 2 V z T D ze -r + V z T D xz V x ) 


This function should be minimized with respect to variable 
V under the constraint, Eq. (5), which is transformed to the 
following condition after discretization 

V T H*d = 0 H = ( S T S ) (23) 


where is matrix with three columns, each corresponding 
to one of the constraints of Eq. (5). The set of columns 
are determined by the kernel (null-space) of matrix E (for 
more details see Ref. 7). This means 

E*d = 0 (24) 


Let us suppose that the set of columns ^ is normalized in 
such a way that 

(25) 


The Euler equation for the problem posed by Eqs. (22) 


and (23) is 



EV+Dh € 'y = H'bdV 


(26) 


where /x is the column matrix of Lagrange multipliers for the 
constraint in Eq. (23). By pre-multiplying Eq. (26) by 
one can prove that 


n = vhDuy 

Now the equation for V, Eq. (26), is rewritten as 


(27) 


in which the following definitions are introduced 

E=([r h sfD[r h s']) 

Du= (M T D[r e ]) 

D hx = (\r h sfD[r t s]) 

A / X ( 21 ) 

Dee = ({r c ] T Dtr t }) 

D z e = ((T/Sf £>[T € ]) 

£>« = ([: r t s] T D[r t s]) 

Classical Considerations 

According to the variational- asymptotical procedure, in 
order to get the next approximation, one should retain only 
the leading energy term with respect to the small parameter 
that contains the unknown functions and the leading inter- 
section term between the unknown function and the rest of 
the functional (for more details see Ref. 7). 

We are left with the following expression 


(j) EV = -(I-H*et1%)Du 7 (28) 

Since E has a kernel, E~ l does not exist. However, the 
pseudo-inverse of E t E^ satisfied the following relations 

eb$-i-hW5 

E^E -I- ^ 

e+ee+ = e+ 

can be foimd (see the Appendix) and the solution of Eq. (28) 
can be represented as 

V=-hE+D h€ l (30) 

Substituting the solution, Eq. (30), into the discretized 
strain energy density, Eq. (20), and keeping only terms with 
the lowest order, which are equal to h° = 1 one obtains 

2t7 = 7 T Ay (31) 


(i)V^v+(i)jWD„ 7 <»> " ,h >1 k D.. - [a.TOa.1 (32) 

The third property, Eq. (29), is taken into account here. 


This is the identical to the classical result for the strain 
energy of laminated plates. 


New Degrees of Freedom 

In order to make our plate functional more flexible with 
respect to the variable x, let us introduce the new unknown 
plate functions such that 

V(x) = %<i(x) + W(x) (33) 


where q is a column matrix of several new unknown func- 
tions, % is a matrix, each column of which represents a 
(-shape form associated with one of the new unknown func- 
tions $(x), which will be named “new degrees of freedom,” 
and W is the new warping to be found 

We will suppose that matrix % is normalized in such 
a way that 

*JH% = / (34) 

The following constraint for W will make the splitting, 
Eq. (33), unique 

W t H% = 0 (35) 


The order of functions g(x) with respect to h may be 
arbitrary *rd it will be supposed to be equal to h° s 1 . 

As a (-shape form of new degrees of freedom let us take 
eigenvectors of matrix E which correspond to the several 
lowest non-zero eigenvalues. Such a matrix % will satisfy 
the following equation 


E%=H%A q (36) 


where is a diagonal matrix of eigenvalues 


[Ai 0 

0 A 2 


L o o 



(37) 


The constraint of Eq. (23), which still might be satisfied 
by TV, ran be combined with the constraint of Eq. (35) after 
introducing the matrix ] hi such a wa Y that 

W t H% = 0 (38) 


where matrix A u is 


•o 

0 

0 

0 

... 0 ' 

0 

0 

0 

0 

... o 

0 

0 

0 

0 

0 

0 

0 

0 

Ai 

... 0 

.0 

0 

0 

0 

A *, - 


( 40 ) 


Calculation of Strain Energy 

Let us assume that we have the correct expansion of V 
through order h 2 

V= V 0 + AVi + /i 2 V 3 (41) 

where Vo denotes the first term of Eq. (33) 

V 0 = %q (42) 

The vector Vq satisfies the equation 

EV 0 = H%A^q (43) 

and vectors Vi and V 2 have to satisfy the constraint found 
in Eq. (38). 

If we have an asymptotically correct expansion for 
Eq. (41), we can calculate an asymptotically correct energy 
for order h° = 1 

2 U= (1) V?EV 0 + 

+ (i) 2 y n r r EVj + A*7 + D hl Vo, x \ + ^ 

+ (l){lf EV 1 + 2 V?EV 0 + 

+ 2 V, T [Du~i + Dhs V 0iI ] + 2 V^ x D zh Vo+ 
+y T D te 7 + 2 V^D zc y + V£ X D XX V 0 , X } 


The underlined terms are equal to zero here because of 
Eq. (43) for V 0 and the constraint of Eq. (38) for V* and V 2 . 
This means we do not need to know the second approxima- 
tion for V (i.e., V 2 ) in order to calculate the energy for order 
h°. 

We shall minimize the functional 
VfEV i + 2 V, r [Du 7 + (D hz - D: h )V 0 , z ] (45) 


Analogously, Eq. (36) can be rewritten as 
E% =H% A u 


in order to find V\ z 

The notation D* h means 

D: h ^[{D hzi ) T ( D hzi f } ( 46 ) 

which comes from 2 V£ D xh V 0 in the fourth line of Eq. (44) 
after integration by parts with respect to x\ and x 2 . 
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The Euler equation for the fir>rt>onal of Eq. (45) is 
EVi + Dhel + {Dhx D xk )\o wX — H% flu (47) 


where fi u is the Lagrange multiplier winch enforces the con- 
straint in Eq. (38). 

Applying a procedure similar to the one used for the 
classical case we can calculate the Lagrange multiplier fiu 

liu = [Dk. 7 + (Dk* - D- zh )V 0rI ] (48) 

and represent the solution of Eq. (47) by 

Vi = Ei [D^ + (Alt - Dlk)Vo,r] (49) 


where the matrix E+ can be found with the following prop- 


erties 


EE% = 


= (50) 


EZEEZ = EZ 


See the Appendix for an explanation of how to calculate 
the matrix 

Substituting this expression into the energy, Eq. (44), 
one obtains 


2 U = V?EV 0 + Q 2 V; T [Dh*7 + D ht V 0 , x ] + 

+ (1){7 T A«7 + 2 1&P«7 + V£ Z P ZZ V 0 ' X } 


where the following notations are Produced 


(51) 


A„ i D tt - [D^fE: [D^] 

P z{ = D zt - [D hz - [D^] ( 52 ) 

Pxx = D zz - [D hz - Lr zh ] T E: [D kl - D' zh ] 


Finally, after substituting the expression for V 0 , 
Eq. (42), the strain energy can be written as 

b | (53) 

where 

A qq = =A, 

A ?I ^ ( *?D hzi % if D^, % ) = A q J 
V = *? D ht = (54) 

a \% T P IlZl % V P *>»A ' 

" *?P*,*,% J 

A tz = [P tIi % P t „%}=A x J 


f 7 V 

A« 


-4« 

b > 

Aqc 

-V 

-V 

l<7.x J 

A zl 

-4x, 

Axx 


Eq. (53) represents the strain energy of a laminated plate 
undergoing deformation which is constrained only in the 
sense that the strain is smalL Displacement and rotation 
of the normal line element appear nonlinearly in the expres- 
sions for 7. On the other hand, the new degrees of freedom 
give rise to simple lm^^r Euler equations. Since the dis- 
placement field is now completely specified in terms of 7, q , 
and q %xy it becomes a simple matter to recover strain and 
stress throughout the plate by use of Eqs. (10) and (17). 
The classical plate energy can be obtained from Eq. (53) 
after it has been minimized with respect to variables q with 
partial derivatives q^ being equal to zero. 


Numerical Results 

Exact solutions for the Stunn-Liouville problem can be 
obtained with symbolic manipulation software, but we were 
only able to carry this out for one- and two-layer plates. 
Thus, we turned to a finite element solution. In Ref. 15 
results for Sturm-Liouville problems with discontinuous co- 
efficients were obtained which agree quite well with the ex- 
act solution. Our finite element method is similar except 
that the orthogonal Jacobi-polynomial-based shape func- 
tions presented in Ref- 14 are used. This means that our 
finite elements have interior degrees of freedom which can 
be added without generation a new element geometry. This 
mesh through the thickness can be as fine as we wish, but 
we must at least have element breaks where discontinuities 
exist. The highest derivatives invoved in the problem are 
first derivatives with respect to and thus continuous 
shape functions can be used. See Ref. 14 for details. 

In this section we present some numerical results for the 
solutions of the Stunn-Liouville problem, compared when 
possible with the exact solution. We start with an error 
analysis of the eigenvalues for the two-layer case, and we 
conclude with the eigenfunctions of one-, two-, and four- 
layer plates. The four-layer results include both a symmetric 
plate and a non-symmetric cross-ply plate. 

For the purpose of discussion only, we align our coordi- 
nate axes with X\ along the length ("longitudinal inplane”) 
and X 2 along the width ("lateral inplane”); these are not 
necessarily aligned with any material axes. For our exam- 
ples, we choose a fiber reinforced composite material which 
has the following material properties 

El = 25 x 10^p6i Et — 10 6 psi 
G L t = 0.5 x lO^psi Gtt — 0.2 x 10 6 psi 
1 *lt — t/ TT — 0.25 

where signifies the direction parallel to the fibers and T the 
transverse direction. These properties, along with a ply an- 
gle, allow the calculation of the matrix D. Note, however, 
that in the example problems and indeed in all laminated 
plates, certain terms in D will vanish. Our theory does 
not require any special terms in D to be zero. Because 
of the vanishing terms in the example problems, there is 
no coupling between out -of- plane and inplane displacement 
components. Thus, certain modes will be entirely inplane 
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and others entirely out -of- plane. Also, since Dqs is constant 
through the thickness for such plates, the thickness mode 
with the lowest eigenvalue for all the examples will be a sine 
function. 

As examples we consider the following four lay-ups: 


[15 c ] 
[- 15715 =] 
[30°/ -30*],^ 
[0790° /0°/90°] 


‘one-Layer” 
•two-layer” 
•four-layer" 
‘cross- ply" 


where the words in quotation marks indicate the terms we 
use to designate the plate under consideration. 


Eigenvalues for Two- Layer Plate 

All of the plates have three zero eigenvalues. These 
correspond with the classical ‘degrees of freedom” of the 
normal line element, the average displacement components 
of that line element. The eigenvalues appear in the above 
equations as well, and the size of these eigenvalues may have 
a bearing on whether the associated degree of freedom is 
an important one. Our present understanding is that the 
smallest eigenvalues are the ones of interest, but this must be 
confirmed by application of the theory with different choices 
for the degrees of freedom. Here we will present a few of the 
smallest non-zero eigenvalues for the two-layer plate and the 
analysis of the error, since the exact solution is available for 
this case. 

In Table 1 the first four nonzero eigenvalues are shown 
for the two-layer plate, both from the exact solution and 
from our finite element approximation with four elements 
(M=4). The order of the shape functions is varied by chang- 
ing J, the number of Jacobi polynomials used to construct 
the higher-order shape functions. The crudest element has 
J= 0, resulting in linear shape functions. Results are shown 
for J~l (quadratic shape functions) and J=2 (cubic shape 
functions). It is seen that all of the finite element results are 
very close to the exact solution. Furthermore, the higher- 
order shape functions greatly improve the accuracy. 


Exact 

Ninncrical Rooks (M = 4) 


1 = 1 

1 = 2 

2.04599617 

115317042 

2.04710246 

204600108 

4.29843523 


4.30275294 

4.29845836 

8.87418665 

10.84944490 

8.93773416 

8.87567158 

10.57173730 

11.12555748 

10.57715155 

10.57176142 


Table 1: Eigenvalues for 2-layer plate 

A more precise analysis of the error for the lowest 
nonzero eigenvalue is presented in Figs. 2 and 3. In Fig. 2 
the relative error versus the size of the eigenvalue problem 
is essentially a straight line on a log-log scale if J, the num- 
ber of Jacobi polynomials in each element, is held constant. 
This means that the size of the matrix increases as the num- 
ber of elements through the thickness is increased. In Fig. 3 


the relative error versus the size of the eigenvalue problem is 
shown for the case when the mesh is held constant; here the 
size of the matrix increases as the order of element shape 
functions is increased. We will present our results below for 
the eigenfunctions based on elements with cubic shape func- 
tions. To analyze plates with a large number of layers, the 
number of elements must be increased. If only the higher- 
order functions were used, such analyses could become ex- 
pensive. However, as seen in these plots, an increase in the 
number of elements, which would be required because of 
discontinuities in the properties between layers, would bring 
the error down without the need of higher-order shape func- 
tions. 



Fig. 2: Relative error for Ai of 2-layer 
plate versus size of matrix E for constant J 



Fig. 3: Relative error for Ai of 2-layer 
plate versus size of matrix E for constant M 


Eigenfunctions for One- Layer Plate 

Eigenfunctions for the smallest four nonzero eigenval- 
ues of the one-layer plate are presented in Figs. 4a - 4d. 
The finite element results, shown as solid lines, were ob- 
tained using 2 elements with cubic shape functions; and the 
exact solution, shown as dashed lines, is indistinguishable 
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from the finite element results in the plots. Fig. 4a shows 
the eigenfunction for a transverse shear mode dominated by 
the lateral inplane displacement, while Fig. 4b shows the 
corresponding one dominated by the longitudinal inplane 
displacement. Note in Fig. 4a the two displacement compo- 
nents are of opposite sign while in Fig. 4b they are of the 
same sign. Fig. 4c shows a higher mode of the same type 
as shown in Fig. 4a. The fourth nonzero eigenvalue has an 
eigenfunction with only out-of-plane displacement. We note 
the smooth character of these functions, winch are dose to 
sines and cosines. Plate theories based on expansion of the 
displacement in power series or trigonometric should provide 
excellent predictive capability for homogeneous plates such 
as this one. 

Eigenfunctions for Two-Layer Plate 

Eigenfunctions for the smallest four nonzero eigenvalues 
of the two-layer plate are presented in Figs. 5a - 5d, along 
with the corresponding eigenvalues. The finite element re- 
sults were obtained using 4 elements with cubic shape func- 
tions, and the exact solution is again indistinguishable from 
the finite element results in the plots. Fig. 5a shows the 
eigenfunction for a transverse shear mode dominated by the 
lateral inplane displacement. Fig. 5b shows the correspond- 
ing mode having more longitudinal inplane displacement, 
but not so strongly dominated by it. Fig. 5c shows & higher 
mode of the same type as shown in Fig. 5a. The fourth 
nonzero eigenvalue again has an eigenfunction with only out- 
of-plane displacement, as shown in Fig. 5d. 

Although the displacement functions shown in Figs. 5c 
and 5d both appear to be smooth, the displacement has a 
discontinuous slope in both Figs. 5a and 5b, in sharp con- 
trast to results for the one-layer plate. Power series approx- 
imations through the entire thickness of the plate would not 
be able to capture this behavior, but the above theory shows 
that the warping induced by this type of displacement is a 
higher-order effect. Thus, a plate theory with these degrees 
of freedom should be a noticeable improvement over classical 
theory. Furthermore, such a theory should be an improve- 
ment over the theory of Refe. 9 and 10. 

Eigenfunctions for Four-Layer Plate 

Eigenfunctions for the smallest four nonzero eigenval- 
ues of the four-layer plate are presented in Figs. 6a - 6d. 
These results were obtained using 8 elements with cubic 
shape functions. Figs. 6a and 6b show the eigenfunctions 
for coupled transverse shear modes, the former having more 
lateral inplane displacement while the latter has more longi- 
tudinal inplane displacement. Fig. 6c shows a higher mode 
of the same type as shown in Fig. 6a. The fourth nonzero 
eigenvalue again has a smooth eigenfunction with only out- 
of-plane displacement, shown in Fig. 6d, as expected. Again, 
in sharp contrast to results for the one-layer plate, the dis- 
placement has discontinuous slopes in Figs. 6a - 6c. The 
symmetry of the plate is exhibited in the symmetry of the 
modes. 





ThiclDtMl 


Fig. 4b: A= 4.9348022 








Fig. 4c. A= 7.89568352 



Fig. 4d: A= 10.57173733 
Fig. 4: Eigenfunctions for 1 -layer plate 
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Fig. 7a: A= 2.78541496 
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Fig. 7b: A= 2.78541496 



Fig. 7c: A= 10.57176142 




Fig. 7d: A= 10.57176142 
7: Eigenfunctions for cross-ply plate 


Eigenfunctions for Crossr-Ply Plate 

Eigenfunctions for the smallest four nonzero eigenval- 
ues of the cross-ply plate are presented in Figs. 7a - 7d. 
These results were obtained using 8 elements with cubic 
shape functions. Figs. 7a and 7b show the eigenfunctions 
for coupled transverse shear modes (with the same eigen- 
values), one having more lateral inplane displacement while 
the other has more longitudinal inplane displacement. The 
third, fourth, and fifth nonzero eigenvalues are also the same. 
The eigenfunctions corresponding to two of these eigenval- 
ues are shown in Figs. 7c and 7d. The out-of-plane mode 
is shown in Fig. Jc 9 having the same smooth eigenfunction 
as before. Fig. 7d shows & higher mode of the same type 
as shown in Fig. 7b. The other eigenfunction for this triple 
root (not shown) is a higher mode of the type as shown 
in Fig. 7a. Again, in contrast to results for the one-layer 
plate, the displacement has discontinuous slopes this time 
in Figs. 7a, 7b, and 7d. The lade of symmetry of the plate 
is exhibited in a lack of symmetry in the biplane modes. 

It should be dear that choosing a priori displacement 
fields which exhibit the character of the inplane modes 
shown in Figs. 4-7 would be virtually impossible. They 
have the appropriate symmetry or asymmetry, as well as the 
discontinuities which reflect the layup. 

Conducting Remarks 

A geometrically nonlinear theory for laminated plates 
is presented based on a combination of the variational- 
asymptotical method and the method of Ritz. The dis- 
placement field is described in terms of the average dis- 
placement of the normal line element and a small number 
of additional functions of the in-plane coordinates of the 
plate. The through-the-thickness shape functions for these 
new “degrees of freedom" are not analytic functions for ar- 
bitrarily laminated plates. Rather, they are eigenfunctions 
of a cer tain Stunn-Liouville problem based on the thickness 
coordinate of the plate. Unlike power series formulations, 
this allows for the correct treatment the known jumps in 
the stress and strain fields. Unlike layerwise variable theo- 
ries, the present theory has only a small number of variables 
in addition to those found in classical plate theory, a num- 
ber which does not depend on the number of layers in the 
plates. Additional equilibrium equations for the plate the- 
ory associated with the new degrees of freedom are simple, 
linear equations - even for a large-displacement theory. 

Since analytical solutions of the Stunn-Liouville prob- 
lem are limited to one- and two-layer plates, an approximate 
finite element solution was obtained. Results obtained for 
these shape functions are presented for a variety of lami- 
nated plate configurations. These results agree well with 
available exact solutions. In the future numerical studies 
will be conducted in order to determine how many and which 
types of degrees of freedom produce the best all-around plate 
theory. 
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